#include "flow.h"
#include <string>
#include <iomanip>
#include <cmath>
// 限制器的需要
double rho_ref, u_ref, p_ref, vol_ref;
// 无需静态， 因为对象不允许复制和移动
size_t flow::imax{0};
size_t flow::ib2{0};

// #define MOPENCL
#define SERIAL
// #define TEST

/**
 * @brief 开辟内存空间
 * 
 * @param imax 包含 ghost 单元在内的单元数
 */
void flow::iniMemory(size_t imax) noexcept
{
    flow::imax = imax;
    flow::ib2 = imax - 1;

    cv.reshape(imax, ncVar);
    cvold.reshape(imax, ncVar);
    p.resize(imax);
    dt.resize(imax);
    diss.reshape(imax, ncVar);
    rhs.reshape(imax, ncVar);
    ls.reshape(imax, ncVar);
    rs.reshape(imax, ncVar);

    buffer1d.resize(imax);
    // cusp 格式有些特殊要求
    buffer2d.reshape(imax + 1, ncVar);

    string opt{inp::options};
    opt += " -D imax=" + std::to_string(imax);
    opt += " -D ib2=" + std::to_string(ib2);
    std::cout << opt << "\n";
    cl_device_type device_type = CL_DEVICE_TYPE_CPU;
    clobj.GetPlatformID(inp::pid);
    clobj.CreateContextFromType(device_type);
    clobj.CreateCommandQueue();
    clobj.CreateProgram(inp::clfile.c_str(), opt.c_str());
    clobj.CreateKernelsInProgram();

    cl_int err;
    cl_context context = clobj.GetContext();
    mem_cv = clCreateBuffer(context, CL_MEM_READ_WRITE,
                            sizeof(double) * imax * ncVar,
                            NULL, &err);
    assert(err == CL_SUCCESS);
    mem_cvold = clCreateBuffer(context, CL_MEM_READ_WRITE,
                               sizeof(double) * imax * ncVar,
                               NULL, &err);
    assert(err == CL_SUCCESS);
    mem_p = clCreateBuffer(context, CL_MEM_READ_WRITE,
                           sizeof(double) * imax,
                           NULL, &err);
    assert(err == CL_SUCCESS);
    mem_dt = clCreateBuffer(context, CL_MEM_READ_WRITE,
                            sizeof(double) * imax,
                            NULL, &err);
    assert(err == CL_SUCCESS);
    mem_diss = clCreateBuffer(context, CL_MEM_READ_WRITE,
                              sizeof(double) * imax * ncVar,
                              NULL, &err);
    assert(err == CL_SUCCESS);
    mem_rhs = clCreateBuffer(context, CL_MEM_READ_WRITE,
                             sizeof(double) * imax * ncVar,
                             NULL, &err);
    assert(err == CL_SUCCESS);
    mem_ls = clCreateBuffer(context, CL_MEM_READ_WRITE,
                            sizeof(double) * imax * ncVar,
                            NULL, &err);
    assert(err == CL_SUCCESS);
    mem_rs = clCreateBuffer(context, CL_MEM_READ_WRITE,
                            sizeof(double) * imax * ncVar,
                            NULL, &err);
    assert(err == CL_SUCCESS);
    mem_buffer1d = clCreateBuffer(context, CL_MEM_READ_WRITE,
                                  sizeof(double) * imax,
                                  NULL, &err);
    assert(err == CL_SUCCESS);
    mem_buffer2d = clCreateBuffer(context, CL_MEM_READ_WRITE,
                                  sizeof(double) * (imax + 1) * ncVar,
                                  NULL, &err);
    assert(err == CL_SUCCESS);

    mem_area = clCreateBuffer(context, CL_MEM_READ_WRITE,
                              sizeof(double) * imax,
                              NULL, &err);
    assert(err == CL_SUCCESS);
    mem_vol = clCreateBuffer(context, CL_MEM_READ_WRITE,
                             sizeof(double) * imax,
                             NULL, &err);
    assert(err == CL_SUCCESS);
}

/**
 * @brief 初始化流场
 *      使用完全气体方程， 条件：入口总压、总温，出口静压
 *      将管内静压全部为出口静压，再利用总温和总压计算
 * @param area 
 */
// void flow::iniFlow(const vec1d<double> &area) noexcept
void flow::iniFlow(const Grid& grd) noexcept
{
    using inp::cpgas;
    using inp::gam1;
    using inp::gamma;
    using inp::gap1;
    using inp::p01;
    using inp::p2;
    using inp::rgas;
    using inp::t01;
    const vec1d<double>& area = grd.area;
    // iniMemory()
    // double gam1{gamma - 1.0}, gap1{gamma + 1.0};
    // double rgas{cpgas * gam1 / gamma};               ///气体常量 R
    double temp = t01 * pow(p2 / p01, gam1 / gamma); ///< 静温
    double rho{p2 / (rgas * temp)};
    double mach = sqrt(2 * (t01 / temp - 1.0) / gam1);
    double as = sqrt(gamma * p2 / rho); ///< 音速
    double u{as * mach};
    double massFlow{rho * u * area(1)};
    double E{(cpgas - rgas) * t01};
    size_t imax = area.size();

#ifdef MOPENCL
    cl_int err;
    cl_kernel initFlow_ker = clobj.GetKernel("InitFlow");

    err = clEnqueueWriteBuffer(clobj.GetCommandQueue(), mem_area, CL_TRUE,
                               0, sizeof(double) * imax, area.data(),
                               0, NULL, NULL);
    err |= clEnqueueWriteBuffer(clobj.GetCommandQueue(), mem_vol, CL_TRUE,
                                0, sizeof(double) * imax, grd.vol.data(),
                                0, NULL, NULL);
    assert(err == CL_SUCCESS);
    // __kernel void InitFlow(global double* cv, global double* p,
    //                     global double* area, const double rho,
    //                     const double massFlow, const double p2,
    // const double E)
    err = clSetKernelArg(initFlow_ker, 0, sizeof(cl_mem), &mem_cv);
    err |= clSetKernelArg(initFlow_ker, 1, sizeof(cl_mem), &mem_p);
    err |= clSetKernelArg(initFlow_ker, 2, sizeof(cl_mem), &mem_area);
    err |= clSetKernelArg(initFlow_ker, 3, sizeof(double), &rho);
    err |= clSetKernelArg(initFlow_ker, 4, sizeof(double), &massFlow);
    err |= clSetKernelArg(initFlow_ker, 5, sizeof(double), &p2);
    err |= clSetKernelArg(initFlow_ker, 6, sizeof(double), &E);
    assert(err == CL_SUCCESS);
    cl_uint work_dim = 1;
    size_t global_size[1]{imax};
    // size_t local_size[1]{} 让实现自己去选
    err = clEnqueueNDRangeKernel(clobj.GetCommandQueue(), initFlow_ker, work_dim,
                                 NULL, global_size, NULL, 0, NULL, NULL);
    assert(err == CL_SUCCESS);
    clFinish(clobj.GetCommandQueue());

    err = clEnqueueReadBuffer(clobj.GetCommandQueue(), mem_cv, CL_TRUE, 0,
                              sizeof(double) * imax * ncVar,
                              cv.data(), 0, NULL, NULL);
    assert(err == CL_SUCCESS);
    err = clEnqueueReadBuffer(clobj.GetCommandQueue(), mem_p, CL_TRUE, 0,
                              sizeof(double) * imax,
                              p.data(), 0, NULL, NULL);

#ifdef TEST
    vec2d<double> clcv = cv;
    vec1d<double> clp = p;
#endif

#endif

#ifdef SERIAL
    for (size_t i = 0; i < imax; ++i)
    {
        cv(i, 0) = rho * area(i);
        cv(i, 1) = massFlow;
        cv(i, 2) = rho * E * area(i);
        p(i) = p2;
        // std::cout << cv(i,1) << "\n";
    }
#endif
//应该定义宏或 inline 函数检查
#ifdef TEST
    if (clcv == cv)
    {
        std::cout << "\n clcv == cv\n";
    }
    else
    {
        std::cerr << "\n clcv != cv\n";
    }
    if (clp == p)
    {
        std::cout << "\n clp == p \n";
    }
    else
    {
        std::cerr << "\n clp != p \n";
    }
#endif
    // exit(EXIT_FAILURE);
    // 限制器参考量？？？ 需要无量纲量
    rho_ref = rho;
    u_ref = u;
    p_ref = p2;
    vol_ref = 1.0;
}
/**
 * @brief 续算读入之前的解
 * 
 */
void flow::rsolut() noexcept
{
}

/**
 * @brief 边界条件
 *      入口亚音速、出口亚音速或超音速： 使用黎曼不变量设置
 * 
 * @param area 单元中心面积
 * @param ib2    
 */
void flow::bcond(const vec1d<double> &area, const size_t ib2) noexcept
{
    using inp::cpgas;
    using inp::gam1;
    using inp::gamma;
    using inp::gap1;
    using inp::p01;
    using inp::p2;
    using inp::rgas;
    using inp::t01;
    using std::cerr;

    // inlet (subsoinc)  268 269 页

    double ud{cv(1, 1) / cv(1, 0)};
    double cs2{gamma * p(1) / cv(1, 0) * area(1)}; ///< 音速平方
    double c02{cs2 + 0.5 * gam1 * ud * ud};
    // outgoing riemann invariant
    double rinv = ud - 2.0 * sqrt(cs2) / gam1;
    double dis = gap1 * c02 / (gam1 * rinv * rinv) - 0.5 * gam1;
    if (dis < 0.0)
    {
        cerr << "warning(bcond0): discriminant < 0 at inflow boundary \t"
             << dis << "\n";
        dis = 1.0e-20;
    }

    double cb = -rinv * gam1 / gap1 * (1.0 + sqrt(dis));
    double cc02 = (cb * cb / c02 < 1.0) ? cb * cb / c02 : 1.0;
    double tb{t01 * cc02};
    double pb = p01 * pow(tb / t01, gamma / gam1);
    double rhob{pb / (rgas * tb)};
    double ub = sqrt(2.0 * cpgas * (t01 - tb));
    // std::cout << std::setprecision(10) << ud << "\t" << cs2 << "\t" << c02 << "\t"
    //         << rinv << "\t" << dis << "\n";
    // std::cout << cb << "\t" << cc02 << "\t" << tb << "\t"
    //         << pb << "\t" << rhob << "\t" << ub << "\n";

    cv(0, 0) = rhob * area(1);
    cv(0, 1) = rhob * area(1) * ub;
    cv(0, 2) = area(1) * (pb / gam1 + 0.5 * rhob * ub * ub);
    p(0) = pb;

    // std::cout << cv(0, 0) << "\t" << cv(0, 1) << "\t"
    //         << cv(0, 2) << "\t" << p(0) << "\n";
    // outlet  subsonic and supersonic

    double rho{cv(ib2 - 1, 0) / area(ib2 - 1)};
    double u{cv(ib2 - 1, 1) / cv(ib2 - 1, 0)};
    double cs = sqrt(gamma * p(ib2 - 1) / rho);

    if (u < cs)
    { // 亚音速流
        pb = p2;
        rhob = rho + (pb - p(ib2 - 1)) / (cs * cs);
        ub = u + (p(ib2 - 1) - pb) / (rho * cs);
    }
    else
    { // 超音速流
        pb = p(ib2 - 1);
        rhob = rho;
        ub = u;
    }
    // imax = ib2+1
    cv(ib2, 0) = rhob * area(ib2);
    cv(ib2, 1) = rhob * area(ib2) * ub;
    cv(ib2, 2) = area(ib2) * (pb / gam1 + 0.5 * rhob * ub * ub);
    p(ib2) = pb;

    // std::cout << cv(ib2, 0) << "\t" << cv(ib2, 1) << "\t"
    //         << cv(ib2, 2) << "\t" << p(ib2) << "\n";
    // exit(EXIT_FAILURE);
}

/**
 * @brief 计算时间步长
 *      谱半径计算有些出人意料
 * 
 * @param vol 
 * @param area 
 * @param x 
 */
void flow::timeStep(const vec1d<double> &vol, const vec1d<double> &area,
                    const vec1d<double> &x) noexcept
{
    using inp::gamma;
    double cs, dx, sprad;
    for (size_t i = 1; i < ib2; ++i)
    {
        cs = sqrt(gamma * p(i) / cv(i, 0) * area(i));
        dx = 0.5 * (x(i + 1) - x(i - 1));
        sprad = cs * sqrt(dx * dx + area(i) * area(i)) + area(i) * std::abs(cv(i, 1) / cv(i, 0));
        dt(i) = vol(i) / sprad;
        // std::cout << std::setprecision(16) << std::scientific
        //         << dt(i) << "\t" << vol(i) << "\n";
        //         << diss(i, 2) << "\n";
    }
    // exit(EXIT_FAILURE);
    dt(0) = dt(1);
    dt(imax - 1) = dt(imax - 2);
}

/**
 * @brief 计算人工粘性项
 *      注意小小的外推和求值方式
 * @param vol 
 * @param beta 
 */
void flow::dissip(const vec1d<double> &vol, double beta) noexcept
{
    using inp::vis2;
    using inp::vis4;
    vec1d<double> &p_sensor = buffer1d;
    vec2d<double> &prediss = buffer2d;
    // vec2d<double>
    for (size_t i = 1; i < ib2; ++i)
    {
        p_sensor(i) = std::abs(p(i + 1) - 2 * p(i) + p(i - 1)) / (p(i + 1) + 2 * p(i) + p(i - 1));
        // std::cout << std::setprecision(16) << std::scientific << p_sensor(i) << "\t"
        // << std::abs(3.5) << "\n";
    }
    // exit(EXIT_FAILURE);
    p_sensor(0) = p_sensor(1);
    p_sensor(imax - 1) = p_sensor(imax - 1);

    // disspation fluxes at i+1/2
    size_t im1, ip1, ip2;
    double eps2, eps4, eval, p_sensor_max;
    // 注意， unsinged 发生了溢出
    for (int i = 0; i < ib2; ++i)
    {
        // std::cout << i << "\n";
        // max
        im1 = (i - 1 < 0) ? 0 : i - 1;
        ip1 = i + 1;
        // min
        ip2 = (i + 2 < (imax - 1)) ? i + 2 : (imax - 1);

        eval = 0.5 * (vol(i) / dt(i) + vol(ip1) / dt(ip1));

        p_sensor_max = (p_sensor(i) < p_sensor(ip1)) ? p_sensor(ip1) : p_sensor(i);

        eps2 = vis2 * p_sensor_max * eval;
        eps4 = vis4 * eval;
        eps4 = (0.0 < eps4 - eps2) ? eps4 - eps2 : 0.0;

        prediss(i, 0) = eps2 * (cv(ip1, 0) - cv(i, 0)) -
                        eps4 * (cv(ip2, 0) - 3 * cv(ip1, 0) +
                                3 * cv(i, 0) - cv(im1, 0));

        prediss(i, 1) = eps2 * (cv(ip1, 1) - cv(i, 1)) -
                        eps4 * (cv(ip2, 1) - 3 * cv(ip1, 1) +
                                3 * cv(i, 1) - cv(im1, 1));

        prediss(i, 2) = eps2 * (cv(ip1, 2) - cv(i, 2)) -
                        eps4 * (cv(ip2, 2) - 3 * cv(ip1, 2) +
                                3 * cv(i, 2) - cv(im1, 2));
    }

    // hybrid multistage schemes
    double beta1{1.0 - beta};
    for (size_t i = 1; i < ib2; ++i)
    {
        diss(i, 0) = beta * (prediss(i, 0) - prediss(i - 1, 0)) + beta1 * diss(i, 0);
        diss(i, 1) = beta * (prediss(i, 1) - prediss(i - 1, 1)) + beta1 * diss(i, 1);
        diss(i, 2) = beta * (prediss(i, 2) - prediss(i - 1, 2)) + beta1 * diss(i, 2);

        // std::cout << std::setprecision(16) << std::scientific
        //         << diss(i, 0) << "\t" << diss(i,1) << "\t"
        //         << diss(i, 2) << "\n";
    }
    // exit(EXIT_FAILURE);
}

/**
 * @brief 计算对流通量和残差
 *      守恒变量的平均
 * 
 * @param area 
 */
void flow::fluxCen(const vec1d<double> &area) noexcept
{
    // 守恒变量的平均（还有压强）， 其它一切物理量均由其推导
    double rhoav, rhouav, rhoEav;
    double rhoHav, pav, areaav, qs;
    vec2d<double> &flux = buffer2d;

    for (size_t i = 0; i < ib2; ++i)
    {
        areaav = 0.5 * (area(i) + area(i + 1));

        rhoav = 0.5 * (cv(i, 0) / area(i) + cv(i + 1, 0) / area(i + 1));
        rhouav = 0.5 * (cv(i, 1) / area(i) + cv(i + 1, 1) / area(i + 1));
        rhoEav = 0.5 * (cv(i, 2) / area(i) + cv(i + 1, 2) / area(i + 1));

        pav = 0.5 * (p(i) + p(i + 1));
        rhoHav = rhoEav + pav;
        qs = rhouav * areaav / rhoav;

        flux(i, 0) = rhoav * qs;
        flux(i, 1) = rhouav * qs + pav * areaav;
        flux(i, 2) = rhoHav * qs;
    }

    // flux + dissipation = rhs
    for (size_t i = 1; i < ib2; ++i)
    {
        rhs(i, 0) = flux(i, 0) - flux(i - 1, 0) - diss(i, 0);
        rhs(i, 1) = flux(i, 1) - flux(i - 1, 1) - diss(i, 1);
        rhs(i, 2) = flux(i, 2) - flux(i - 1, 2) - diss(i, 2);
        // std::cout << std::setprecision(16) << std::scientific
        //         << rhs(i, 0) << "\t" << rhs(i,1) << "\t"
        //         << rhs(i, 2) << "\n";
    }
    // exit(EXIT_FAILURE);
}

void flow::SrcTerm(const vec1d<double> &area) noexcept
{
    double da;
    for (size_t i = 1; i < ib2; ++i)
    {
        da = 0.5 * (area(i + 1) - area(i - 1));
        rhs(i, 1) -= p(i) * da;
        // std::cout << std::setprecision(16) << std::scientific
        //         << rhs(i, 1) << "\n";
        //         << diss(i, 2) << "\n";
    }
    // exit(EXIT_FAILURE);
}

/**
 * @brief 残值光顺
 * 
 */
void flow::Irsmoo() noexcept
{
    using inp::epsirs;

    vec1d<double> &d = buffer1d;
    d(0) = 0.0;
    rhs(0, 0) = 0.0;
    rhs(0, 1) = 0.0;
    rhs(0, 2) = 0.0;

    double t, eps2{2.0 * epsirs + 1.0};
    for (size_t i = 1; i < ib2; ++i)
    {
        t = 1.0 / (eps2 - epsirs * d(i - 1));
        d(i) = t * epsirs;
        rhs(i, 0) = t * (rhs(i, 0) + epsirs * rhs(i - 1, 0));
        rhs(i, 1) = t * (rhs(i, 1) + epsirs * rhs(i - 1, 1));
        rhs(i, 2) = t * (rhs(i, 2) + epsirs * rhs(i - 1, 2));
        // std::cout << std::setprecision(16) << d(i) << "\n";
        // std::cout << std::setprecision(16) << rhs(i, 0) << "\n";
    }

    size_t i = ib2 - 1;
    for (size_t j = 3; j <= ib2; ++j)
    {
        --i;
        rhs(i, 0) = rhs(i, 0) + d(i) * rhs(i + 1, 0);
        rhs(i, 1) = rhs(i, 1) + d(i) * rhs(i + 1, 1);
        rhs(i, 2) = rhs(i, 2) + d(i) * rhs(i + 1, 2);
        // std::cout << std::setprecision(16) << rhs(i, 0) << "\n";
    }
}

void flow::UpdateConservativeVar(const vec1d<double> &area) noexcept
{
    using inp::gam1;

    for (size_t i = 1; i < ib2; ++i)
    {
        cv(i, 0) = cvold(i, 0) - rhs(i, 0);
        cv(i, 1) = cvold(i, 1) - rhs(i, 1);
        cv(i, 2) = cvold(i, 2) - rhs(i, 2);

        p(i) = gam1 * (cv(i, 2) / area(i) - 0.5 * cv(i, 0) / area(i) * (cv(i, 1) / cv(i, 0) * cv(i, 1) / cv(i, 0)));
    }
}

void flow::Conver(double &drho, const string &convergence_file_name, const vec1d<double> &area) noexcept
{
    using inp::gamma;

    drho = 0.0;
    double drmax{0.0}, av_mass_flow{0.0};
    double dr{0.0}, rho{0.0}, u{0.0}, c{0.0};
    size_t index_drmax{0}, nsup{0};

    for (size_t i = 1; i < ib2; ++i)
    {
        dr = cv(i, 0) - cvold(i, 0);
        av_mass_flow += cv(i, 1);
        drho += dr * dr;
        if (std::abs(dr) > drmax)
        {
            drmax = std::abs(dr);
            index_drmax = i;
        }

        rho = cv(i, 0) / area(i);
        u = cv(i, 1) / cv(i, 0);
        c = sqrt(gamma * p(i) / rho);
        if (c < u)
            ++nsup;
    }
    av_mass_flow /= (ib2 - 1);

    static double drho1{0.0};
    static size_t iter{0};
    ++iter;
    // 以第一步平均密度差值作为基准
    if (iter == 1)
        drho1 = sqrt(drho / (ib2 - 1));

    drho = sqrt(drho / (ib2 - 1)) / drho1;

    std::cout << std::scientific << std::setprecision(9)
              << iter << "\t" << log10(drho) << "\t"
              << drmax << "\t" << index_drmax << "\t"
              << av_mass_flow << "\t" << cv(1, 1) - cv(ib2 - 1, 1) << "\t"
              << nsup << "\n";
    // 不应该这样干
    // using inp::convergence_file;
    // std::ofstream convergence_file{inp::fnconv, std::ios_base::out};
    // convergence_file << std::scientific << std::setprecision(9)
    //                  << iter << "\t" << log10(drho) << "\t"
    //                  << drmax << "\t" << index_drmax << "\t"
    //                  << av_mass_flow << "\t" << cv(1, 1) - cv(ib2 - 1, 1) << "\t"
    //                  << nsup << "\n";
    // // 不应该这样做， 先这样
    // convergence_file.close();
}

void flow::solver(const Grid &grid, double &drho) noexcept
{
    // using namespace inp;
    using inp::ark;
    using inp::betrk;
    using inp::cfl;
    using inp::epsirs;
    using inp::kdissip;
    using inp::ldiss;
    using inp::lsmoo;
    using inp::nrk;

    // store previous solution, set dissipation = 0;
    cvold.copyFrom(cv);
    diss = 0.0;
    // calculate time step
    timeStep(grid.vol, grid.area, grid.x);
    // std::cout << "4\n";
    double fac;
    const vec1d<double> &vol = grid.vol;

    for (size_t i = 0; i != nrk; ++i)
    {
        // std::cout << cv(1, 0) << "\t" << rhs(1, 0) << "\n";
        if (kdissip == "c" || kdissip == "C")
        {
            // do something
            // std::cout << "1\n";

            if (ldiss[i])
                dissip(vol, betrk[i]);
            // std::cout << "6\n";
            // std::cout << cv(1, 0) << "\t" << rhs(1, 0) << "\n";
            fluxCen(grid.area);
            // std::cout << "5\n";
            // std::cout << cv(1, 0) << "\t" << rhs(1, 0) << "\n";
        }
        else if (kdissip == "r" || kdissip == "R")
        {
            LRStateRoe(grid.area, grid.vol);
            FluxRoe(grid.area);
        }
        else
        {
            std::cerr << "now, only support central scheme and roe scheme\n";
            exit(EXIT_FAILURE);
        }

        SrcTerm(grid.area);
        // std::cout << cv(1, 0) << "\t" << rhs(1, 0) << "\n";

        // 每一阶 残差前的系数
        fac = ark[i] * cfl;
        double adtv;
        // std::cout << fac << "\n";
        for (size_t i = 1; i < ib2; ++i)
        {
            adtv = fac * (dt(i) / vol(i));
            rhs(i, 0) *= adtv;
            rhs(i, 1) *= adtv;
            rhs(i, 2) *= adtv;
            // std::cout << std::setprecision(15) << std::scientific << i << "\t"
            //           << rhs(i, 0) << "\n";
        }
        // std::exit(EXIT_FAILURE);
        // std::cout << "sss\n";

        // 残差光顺
        if (lsmoo[i] && (epsirs > 0.0))
            Irsmoo();
        // std::cout << cv(1, 0) << "\t" << rhs(1, 0) << "\n";
        // exit(EXIT_FAILURE);
        // update
        UpdateConservativeVar(grid.area);

        bcond(grid.area, ib2);
    }
        // exit(EXIT_FAILURE);

    Conver(drho, inp::fnconv, grid.area);
}

void flow::LRStateCusp(const vec1d<double> &area) noexcept
{
    vec2d<double> &du = buffer2d;
    // 注意， 移动了一位
    for (int i = 1; i < imax; ++i)
    {
        du(i, 0) = cv(i, 0) / area(i) - cv(i - 1, 0) / area(i - 1);
        du(i, 1) = cv(i, 1) / cv(i, 0) - cv(i - 1, 1) / cv(i - 1, 0);
        du(i, 2) = p(i) - p(i - 1);
    }
    du(0, 0) = du(1, 0);
    du(0, 1) = du(1, 1);
    du(0, 2) = du(1, 2);
    // buffer2d imax+1, 3
    du(imax, 0) = du(ib2, 0);
    du(imax, 1) = du(ib2, 1);
    du(imax, 2) = du(ib2, 2);

    auto CuspLimiter = [](double d1, double d2) {
        return 1.0 - (std::pow((d1 - d2) / (std::abs(d1) + std::abs(d2) + 1.0e-20), 2.0));
    };

    for (int i = 0; i < imax; ++i)
    {
        ls(i, 0) = cv(i, 0) / area(i) + 0.25 * (du(i, 0) + du(i + 2, 0)) * CuspLimiter(du(i, 0), du(i + 2, 0));
        ls(i, 1) = cv(i, 1) / cv(i, 0) + 0.25 * (du(i, 1) + du(i + 2, 1)) * CuspLimiter(du(i, 1), du(i + 2, 1));
        ls(i, 2) = p(i) + 0.25 * (du(i, 2) + du(i + 2, 2)) * CuspLimiter(du(i, 2), du(i + 2, 2));

        rs(i, 0) = cv(i + 1, 0) / area(i + 1) - 0.25 * (du(i, 0) + du(i + 2, 0)) * CuspLimiter(du(i, 0), du(i + 2, 0));
        rs(i, 1) = cv(i + 1, 1) / cv(i + 1, 0) - 0.25 * (du(i, 1) + du(i + 2, 1)) * CuspLimiter(du(i, 1), du(i + 2, 1));
        rs(i, 2) = p(i + 1) - 0.25 * (du(i, 2) + du(i + 2, 2)) * CuspLimiter(du(i, 2), du(i + 2, 2));
    }
}

void flow::FluxCusp(const vec1d<double> &area) noexcept
{
}

// inline double MUSCL3(double af, double bf, double eps)
// {
//     return (bf * (2.0 * af * af + eps) + af * (bf * bf + 2.0 * eps)) /
//             (2.0 * af * af + 2.0 * bf * bf - af * bf + 3.0 * eps + 1.0e-30);
// }

void flow::LRStateRoe(const vec1d<double> &area, const vec1d<double> &vol) noexcept
{
    using inp::iorder;
    using inp::limfac;
    // using inp::p_ref;
    // using inp::rho_ref;
    // using inp::u_ref;
    // using inp::vol_ref;
    vec2d<double> &du = buffer2d;

    // 限制器中的无量纲所要求的
    double limfac3 = limfac * limfac * limfac;
    double rvol_ref = 1.0 / pow(vol_ref, 1.5);
    double eps2rho = limfac3 * rho_ref * rho_ref * rvol_ref;
    double eps2u = limfac3 * u_ref * u_ref * rvol_ref;
    double eps2p = limfac3 * p_ref * p_ref * rvol_ref;
    double deltr[3], deltl[3];
// /**
#ifdef MOPENCL
    cl_int err;
    //     kernel void
    //     LrStateRoe(global const double*cv, global const double*area,
    //             global double* du, global double* rs,
    //             global double* ls, global const double* p,
    //             global const double *vol,
    //             const int iorder, const double eps2rho,
    //             const double eps2u, const double eps2p)
    cl_kernel LrStateRoe_ker = clobj.GetKernel("LrStateRoe");
    err = clEnqueueWriteBuffer(clobj.GetCommandQueue(), mem_cv, CL_TRUE,
                               0, sizeof(double) * imax * ncVar, cv.data(),
                               0, NULL, NULL);
    err |= clEnqueueWriteBuffer(clobj.GetCommandQueue(), mem_p, CL_TRUE,
                                0, sizeof(double) * imax, p.data(),
                                0, NULL, NULL);
    assert(err == CL_SUCCESS);

    err = clSetKernelArg(LrStateRoe_ker, 0, sizeof(cl_mem), &mem_cv);
    err |= clSetKernelArg(LrStateRoe_ker, 1, sizeof(cl_mem), &mem_area);
    err |= clSetKernelArg(LrStateRoe_ker, 2, sizeof(cl_mem), &mem_buffer2d);
    err |= clSetKernelArg(LrStateRoe_ker, 3, sizeof(cl_mem), &mem_rs);
    err |= clSetKernelArg(LrStateRoe_ker, 4, sizeof(cl_mem), &mem_ls);
    err |= clSetKernelArg(LrStateRoe_ker, 5, sizeof(cl_mem), &mem_p);
    err |= clSetKernelArg(LrStateRoe_ker, 6, sizeof(cl_mem), &mem_vol);
    err |= clSetKernelArg(LrStateRoe_ker, 7, sizeof(int), &iorder);
    err |= clSetKernelArg(LrStateRoe_ker, 8, sizeof(double), &eps2rho);
    err |= clSetKernelArg(LrStateRoe_ker, 9, sizeof(double), &eps2u);
    err |= clSetKernelArg(LrStateRoe_ker, 10, sizeof(double), &eps2p);
    assert(err == CL_SUCCESS);

    cl_uint work_dim = 1;
    size_t global_size[1]{ib2};
    // 这个核函数 组内数量不能任选 要么经过精确的分割 要么全给  数据有依赖关系
    // 必须全部塞进去  栅栏只对工作组内的工作项有作用， 工作组之间没法同步
    // 两端数据没法保持有效性，但可以将任务分割，写多个核函数，本kernel对du写了又读了
    size_t local_size[1]{ib2};
    err = clEnqueueNDRangeKernel(clobj.GetCommandQueue(), LrStateRoe_ker,
                                 work_dim, NULL, global_size, local_size, 
                                 0, NULL, NULL);
    assert(err == CL_SUCCESS);
    
    clFinish(clobj.GetCommandQueue());

#ifdef TEST
    // 把下面数据从全局内存中拿出来纯粹是为了测试
    err = clEnqueueReadBuffer(clobj.GetCommandQueue(), mem_ls, CL_TRUE,
                               0, sizeof(double) * imax * ncVar, ls.data(),
                               0, NULL, NULL);
    assert(err == CL_SUCCESS);

    err = clEnqueueReadBuffer(clobj.GetCommandQueue(), mem_rs, CL_TRUE,
                               0, sizeof(double) * imax * ncVar, rs.data(),
                               0, NULL, NULL);
    assert(err == CL_SUCCESS);
    // 注意 buffer2d 的大小
    err = clEnqueueReadBuffer(clobj.GetCommandQueue(), mem_buffer2d, CL_TRUE,
                               0, sizeof(double) * (imax+1) * ncVar, du.data(),
                               0, NULL, NULL);
    assert(err == CL_SUCCESS);

    vec2d<double> clls = ls;
    vec2d<double> clrs = rs;
    vec2d<double> cldu = du;
#endif
#endif
// **/
#ifdef SERIAL
    if (2 <= iorder)
    {
        // 注意错了一位 第一个代表第0个
        for (int i = 1; i < imax; ++i)
        {
            du(i, 0) = cv(i, 0) / area(i) - cv(i - 1, 0) / area(i - 1);
            du(i, 1) = cv(i, 1) / cv(i, 0) - cv(i - 1, 1) / cv(i - 1, 0);
            du(i, 2) = p(i) - p(i - 1);
        }
        du(0, 0) = du(1, 0);
        du(0, 1) = du(1, 1);
        du(0, 2) = du(1, 2);
        du(imax, 0) = du(ib2, 0);
        du(imax, 1) = du(ib2, 1);
        du(imax, 2) = du(ib2, 2);
    }
    // 4.117
    auto MUSCL3 = [](double af, double bf, double eps) {
        return (bf * (2.0 * af * af + eps) + af * (bf * bf + 2.0 * eps)) /
                (2.0 * af * af + 2.0 * bf * bf - af * bf + 3.0 * eps);
            // (2.0 * af * af + 2.0 * bf * bf - af * bf + 3.0 * eps + 1.0e-30);
    };

    double volt{0.0};
    if (iorder == 3)
    {
        for (int i = 0; i < ib2; ++i)
        {
            volt = pow(0.5 * (vol(i + 1) + vol(i)), 1.5);
            deltr[0] = MUSCL3(du(i + 2, 0), du(i + 1, 0), eps2rho * volt);
            deltl[0] = MUSCL3(du(i + 1, 0), du(i, 0), eps2rho * volt);

            deltr[1] = MUSCL3(du(i + 2, 1), du(i + 1, 1), eps2u * volt);
            deltl[1] = MUSCL3(du(i + 1, 1), du(i, 1), eps2u * volt);

            deltr[2] = MUSCL3(du(i + 2, 2), du(i + 1, 2), eps2p * volt);
            deltl[2] = MUSCL3(du(i + 1, 2), du(i, 2), eps2p * volt);

            rs(i, 0) = cv(i + 1, 0) / area(i + 1) - 0.5 * deltr[0];
            rs(i, 1) = cv(i + 1, 1) / cv(i + 1, 0) - 0.5 * deltr[1];
            rs(i, 2) = p(i + 1) - 0.5 * deltr[2];

            ls(i, 0) = cv(i, 0) / area(i) + 0.5 * deltl[0];
            ls(i, 1) = cv(i, 1) / cv(i, 0) + 0.5 * deltl[1];
            ls(i, 2) = p(i) + 0.5 * deltl[2];
        }
    }

    // ls(ib2, 0) = 0.0;
    // ls(ib2, 1) = 0.0;
    // ls(ib2, 2) = 0.0;
    // rs(ib2, 0) = 0.0;
    // rs(ib2, 1) = 0.0;
    // rs(ib2, 2) = 0.0;
#endif

    if (iorder != 3)
    {
        std::cerr << "roe scheme now only suport 3order\n";
        exit(EXIT_FAILURE);
    }



// /**
#ifdef TEST
    if (clls == ls)
    {
        std::cout << "\n clls == ls\n";
    }
    else
    {
        std::cout << "\n clls != ls\n";
        std::cerr << "\n clls != ls\n";
        for(size_t i = 0; i < ls.dim1(); ++i)
        {
            for(size_t j = 0; j < ls.dim2(); ++j)
            {
                if(abs(clls(i, j) - ls(i, j)) > 2e-16 )
                    std::cout << i << "\t" << j << "\t"
                            << clls(i, j) - ls(i, j) << "\n";
            }
        }
    }
    if (clrs == rs)
    {
        std::cout << "\n clrs == rs \n";
    }
    else
    {
        std::cout << "\n clrs != rs\n";
        std::cerr << "\n clrs != rs \n";
        for(size_t i = 0; i < ls.dim1(); ++i)
        {
            for(size_t j = 0; j < ls.dim2(); ++j)
            {
                if(abs(clrs(i, j) - rs(i, j)) > 2e-16 )
                    std::cout << i << "\t" << j << "\t"
                            << clrs(i, j) - rs(i, j) << "\n";
            }
        }
    }

    if (cldu == du)
    {
        std::cout << "\n cldu == du \n";
    }
    else
    {
        std::cout << "\n cldu != du\n";
        std::cerr << "\n cldu != du\n";
        for(size_t i = 0; i < du.dim1(); ++i)
        {
            for(size_t j = 0; j < du.dim2(); ++j)
            {
                // if(clls(i, j) != ls(i, j))
                    // std::cout << i << "\t" << j << "\t"
                            // << cldu(i, j) << "\t" << du(i, j) << "\n";
                std::cout << std::setprecision(16) << std::scientific << cldu(i, j) - du(i, j) << '\t';
            }
            std::cout << "\n";
        }
    }
#endif
// **/
}

// 串行和并行没有出现不相同的情况， 完全相同
void flow::FluxRoe(const vec1d<double> &area) noexcept
{
    using inp::epsentr;
    using inp::gam1;
    using inp::gamma;
    double avs{0.0};
    double avflux[3], lrho, lu, lp, lq, lh, rrho, ru, rp, rq, rh;
    double rhot, ut, ht, ct, q2t, c2t, tmp, fdiss[3];
    double lambda[3], delta, dv, eps1, eps2, eps3;
    vec2d<double> &flux = buffer2d;
/**
#ifdef MOPENCL
    cl_int err;
    // kernel void
    // FluxRoeKer(global const double* area, global const double* ls,
    //         global const double* rs, global double* flux,
    //         global double* rhs, const double gamma,
    //         const double epsentr)
    cl_kernel flux_roe_ker = clobj.GetKernel("FluxRoeKer");
    cl_command_queue command_queue = clobj.GetCommandQueue();
#ifdef TEST
    err = clEnqueueWriteBuffer(clobj.GetCommandQueue(), mem_ls, CL_TRUE,
                               0, sizeof(double) * imax * ncVar, ls.data(),
                               0, NULL, NULL);
    err |= clEnqueueWriteBuffer(clobj.GetCommandQueue(), mem_rs, CL_TRUE,
                                0, sizeof(double) * imax * ncVar, rs.data(),
                                0, NULL, NULL);
    assert(err == CL_SUCCESS);
#endif
    // 不需要向内存对象传递数据
    err = clSetKernelArg(flux_roe_ker, 0, sizeof(cl_mem), &mem_area);
    err |= clSetKernelArg(flux_roe_ker, 1, sizeof(cl_mem), &mem_ls);
    err |= clSetKernelArg(flux_roe_ker, 2, sizeof(cl_mem), &mem_rs);
    err |= clSetKernelArg(flux_roe_ker, 3, sizeof(cl_mem), &mem_buffer2d);
    err |= clSetKernelArg(flux_roe_ker, 4, sizeof(cl_mem), &mem_rhs);
    err |= clSetKernelArg(flux_roe_ker, 5, sizeof(double), &gamma);
    err |= clSetKernelArg(flux_roe_ker, 6, sizeof(double), &epsentr);
    assert(err == CL_SUCCESS);

    cl_uint work_dim = 1;
    size_t global_size[1]{ib2};
    // 这个核函数 组内数量不能任选 要么经过精确的分割 要么全给  数据有依赖关系
    size_t local_size[1]{ib2};
    err = clEnqueueNDRangeKernel(command_queue, flux_roe_ker, 1, NULL,
                                global_size, local_size, 0, NULL, NULL);
    assert(err == CL_SUCCESS);
    clFinish(command_queue);
    // 没有把全部内存全取出来， 最后面的用不到
    err = clEnqueueReadBuffer(command_queue, mem_rhs, CL_TRUE, 0, 
                        sizeof(double)*ncVar*ib2, rhs.data(), 0,
                        NULL, NULL);
    assert(err == CL_SUCCESS);

#ifdef TEST
    err = clEnqueueReadBuffer(command_queue, mem_buffer2d, CL_TRUE, 0, 
                        sizeof(double)*ncVar*ib2, flux.data(), 0,
                        NULL, NULL);
    assert(err == CL_SUCCESS);

    vec2d<double>& clflux = flux;
    vec2d<double>& clrhs = rhs;
#endif

#endif
**/


#ifdef SERIAL
    auto EntropyCorrection = [](double lam, double dlt) {
        return (dlt < lam) ? lam : (lam * lam + dlt * dlt) / dlt * 0.5;
    };

    for (int i = 0; i < ib2; ++i)
    {
        // 面积为单元中心的面积，需要求个平均
        avs = 0.5 * (area(i) + area(i + 1));

        lrho = ls(i, 0);
        lu = ls(i, 1);
        lp = ls(i, 2);
        lq = lu * lrho;
        // 总颔
        lh = gamma / gam1 * lp / lrho + 0.5 * lu * lu;

        rrho = rs(i, 0);
        ru = rs(i, 1);
        rp = rs(i, 2);
        rq = ru * rrho;
        rh = gamma / gam1 * rp / rrho + 0.5 * ru * ru;
        // average of left and right convective fluxes
        avflux[0] = lq + rq;
        avflux[1] = lq * lu + lp + rq * ru + rp;
        avflux[2] = lq * lh + rq * rh;

        // dissipative fluxes
        rhot = sqrt(lrho * rrho);
        tmp = rhot / lrho;
        ut = (lu + ru * tmp) / (1.0 + tmp);
        ht = (lh + rh * tmp) / (1.0 + tmp);
        q2t = ut * ut;
        c2t = gam1 * (ht - 0.5 * q2t);
        ct = sqrt(c2t);

        lambda[0] = abs(ut - ct);
        lambda[1] = abs(ut);
        lambda[2] = abs(ut + ct);
        delta = epsentr * ct;
        lambda[0] = EntropyCorrection(lambda[0], delta);
        lambda[1] = EntropyCorrection(lambda[1], delta);
        lambda[2] = EntropyCorrection(lambda[2], delta);
        dv = ru - lu;
        eps1 = lambda[0] * (rp - lp - rhot * ct * dv) / c2t * 0.5;
        eps2 = lambda[1] * (rrho - lrho - (rp - lp) / c2t);
        eps3 = lambda[2] * (rp - lp + rhot * ct * dv) / c2t * 0.5;

        fdiss[0] = eps1 + eps2 + eps3;
        fdiss[1] = eps1 * (ut - ct) + eps2 * ut + eps3 * (ut + ct);
        fdiss[2] = eps1 * (ht - ct * ut) + eps2 * q2t * 0.5 + eps3 * (ht + ct * ut);
        // 面总的flux
        flux(i, 0) = 0.5 * avs * (avflux[0] - fdiss[0]);
        flux(i, 1) = 0.5 * avs * (avflux[1] - fdiss[1]);
        flux(i, 2) = 0.5 * avs * (avflux[2] - fdiss[2]);
    }

    // 收集通量得到残差
    for (int i = 1; i < ib2; ++i)
    {
        rhs(i, 0) = flux(i, 0) - flux(i - 1, 0);
        rhs(i, 1) = flux(i, 1) - flux(i - 1, 1);
        rhs(i, 2) = flux(i, 2) - flux(i - 1, 2);
    }
#endif
/**
#ifdef TEST
    if (clflux == flux)
    {
        std::cout << "\n clflux == flux\n";
    }
    else
    {
        std::cout << "\n clflux != flux\n";
        std::cerr << "\n clflux != flux\n";
        for(size_t i = 0; i < flux.dim1(); ++i)
        {
            for(size_t j = 0; j < flux.dim2(); ++j)
            {
                if(abs(clflux(i, j) - flux(i, j)) > 2e-16 )
                    std::cout << i << "\t" << j << "\t"
                            << clflux(i, j) - flux(i, j) << "\n";
            }
        }
    }
    if (clrhs == rhs)
    {
        std::cout << "\n clrhs == rhs \n";
    }
    else
    {
        std::cout << "\n clrhs != rhs\n";
        std::cerr << "\n clrhs != rhs \n";
        for(size_t i = 0; i < rhs.dim1(); ++i)
        {
            for(size_t j = 0; j < rhs.dim2(); ++j)
            {
                if(abs(clrhs(i, j) - rhs(i, j)) > 2e-16 )
                    std::cout << i << "\t" << j << "\t"
                            << clrhs(i, j) - rhs(i, j) << "\n";
            }
        }
    }
#endif
**/

}